/*
 *--------------------------------------------------------------
 *	RELAX.C
 *--------------------------------------------------------------
 */

/*
 *-------------------------------------------------------------
 *	 INCLUDE
 *-------------------------------------------------------------
 */	
#ifdef GRAPHICS
   #include <gl/gl.h>
#endif

#include <stdio.h>
#include <math.h>
#include "common.h"


/*
 *-----------------------------------------------------------------
 * 	EXTERN GLOBAL VARIABLES
 *-----------------------------------------------------------------
 */
/* declared in list.c */
extern CELL *headC;
extern CELL *zedC;

/* declared in bucket.c */
extern BUCKET *head[BUCKET_SIZE_X][BUCKET_SIZE_Y];
extern BUCKET *zed[BUCKET_SIZE_X][BUCKET_SIZE_Y];

/* declared in main.c */
extern double 	radius, nradius, rsquare;
extern float 	min_x, max_x, min_y, max_y;
extern int      bx, by, npoints;
extern float 	wd;
extern float 	adhes[HOW_MANY_TYPES][HOW_MANY_TYPES];
extern float 	wa;
extern float 	wrp;
extern float 	sizeX, sizeY;
extern double   AniX, AniY, AniCommon;

extern voronoiType	*myVoronoiList;	/* declared in main.c */
extern faceType		*myFacesList;

extern float	orientation;
extern double 	Ox, Oy, OxSq, OySq;
/*
 *-----------------------------------------------------------------
 * 	LOCAL GLOBAL VARIABLES
 *-----------------------------------------------------------------
 */
double minForce, maxForce;

/*
 *-----------------------------------------------------------------------
 * This function goes over all points present in the
 * buckets and computes a new position for them according
 * to a relaxation procedure.
 *-----------------------------------------------------------------------
 */
void relaxation(int steps)
{
  CELL *cell;
  int i, ix, iy;
  double dxi, dyi, f;


#ifdef GRAPHICS
  display(FALSE);
#endif

  minForce = 1e20, maxForce = -1e20;

  for(i = 0; i < steps; i++)
  {
    /* Compute all forces for each cell */
    cell = headC->next;
    while( cell != zedC )
    { /* while there is some element in the list */
      dxi = dyi = 0;
      compute_forces( cell, &dxi, &dyi );
      cell->fx = dxi;
      cell->fy = dyi;
      // f = sqrt(fx*fx + fy*fy); 
      // if( f < minForce ) minForce = f; 
      // if( f > maxForce ) maxForce = f; 

      /* go to the next cell in the global list */
      cell = cell->next;
    }
    if (steps > 50){
      if (i % 50 == 0)
	  printf("step = %d\n",i);
    }
    /*fprintf( stderr, "min disp = %lg max disp = %lg\n",
	     minForce, maxForce );*/
    /* Update the new positions for the cells */
    cell = headC->next;
    while( cell != zedC ) { /* while there is some element in the list */
      
      /* find the current bucket for this cell */
      findBucket( cell->x, cell->y, &ix, &iy );
      
      /* finally update the new position */
      cell->x += cell->fx;
      cell->y += cell->fy;
      
      /* check if new position is out of boundaries */
      fixCoordinates( &(cell->x), &(cell->y) );
      
      checkBuckets( cell, ix, iy );
      
      /* go to the next cell in the global list */
      cell = cell->next;
    }
    
#ifdef GRAPHICS
    display(FALSE);
#endif

  } /* end for number of relaxation steps */
}

/*
 *-------------------------------------------------------------
 * This function checks if a cell has changed its bucket and
 * if so,  update the bucket structure
 *-------------------------------------------------------------
 */
void checkBuckets(CELL *s,  int cbx,  int cby)
{
  int i, j;
  BUCKET *t, *tr;

  /* find the new bucket */
  findBucket(s->x, s->y, &i, &j);
    
  if( (i > bx) || (j > by) || (i < 0 ) || (j < 0))
      errorMsg("Trying to reach outside a bucket...");

  if((cbx != i) || (cby != j)){ /* mean that cell has changed buckets */
    t = balloc();
    tr = (head[cbx][cby])->next;
    while(tr->info != s) tr = tr->next;
    t->info = tr->info;
    insert_bucket(t,head[i][j],zed[i][j]);
    remove_bucket(tr,head[cbx][cby],zed[cbx][cby]);
  }
}


/*
 *-----------------------------------------------------------------------
 * This function receives a pointer to the cell to relax. The relaxation
 * procedure is based on the relaxation radius plus the current
 * cell position. The displacement in x and y is returned
 * through sumdx and sumdy. The neighbors are all the cells which
 * fall inside the given repulsive radius
 *-----------------------------------------------------------------------
 */
void compute_forces(CELL *cell, double *sumdx, double *sumdy)
{
  BUCKET *bucket;
  int i, j, ii, jj, b_cell_x, b_cell_y;
  double dxi = 0.0, dyi = 0.0;
  double x, y;

  *sumdx = *sumdy = 0.0;

  /* find the current bucket for this cell */
  findBucket( cell->x, cell->y, &b_cell_x, &b_cell_y );

  for(ii = b_cell_x-1; ii <= b_cell_x+1; ii++)
  {
    for(jj = b_cell_y-1; jj <= b_cell_y+1; jj++)
    {
      /* we want toroidal boundary conditions */
      i = (ii + bx) % bx;
      j = (jj + by) % by;
      
      /* get the first cell on this bucket */
      bucket = (head[i][j])->next;

      while( bucket != zed[i][j] )
      { /* go over all cells in this bucket */
        /* check if the cell in the bucket is not the cell itself */
        if ( cell != (bucket->info) )
        {
            x = (bucket->info)->x;
            y = (bucket->info)->y;
	  
            /* this part has to do with the cells at the borders */
            if ( ii != i ) {
                if( ii > 0 ) x += max_x;
                else x -= max_x;
            }
            if( jj != j ) {
                if ( jj > 0 ) y += max_y;
                else y -= max_y;
            }
            /* get the force between cells */
            force( cell, (bucket->info)->ctype, x, y, &dxi, &dyi );
            *sumdx += dxi;
            *sumdy += dyi;
        }

        /* go to the next cell in this bucket */
        bucket = bucket->next;
      }
    }
  }
}



/*
 *---------------------------------------------------------------------
 * Compute the force between two cells. Returns the x and y
 * displacement components through the dix and diy variables
 * Returns TRUE if the computation is valid and FALSE if not
 * Not valid computation means that
 * 	- the distance between the two cells is more than the
 *	  repulsive radius distance
 *---------------------------------------------------------------------
 */
int force(CELL *s, CELL_TYPE ttype, double xt, double yt, double *dxi, double *dyi)
{
  double d, dsquare, f;
  double dx, dy, adh;
  double x, y, csquare, g=1.0;
  double tmpx, tmpy;

  *dxi = 0.0;
  *dyi = 0.0;

  dx = s->x - xt;
  dy = s->y - yt;
	
  dsquare = dx*dx + dy*dy;
	
  /* This part is a test to see if we could simulate the
     effect of an "area of different influence" in the
     domain. In this case the area is a circle which
     affect the value of "g" used later. For now this
     feature is disabled and "g" remains constant and
     equal to 1.0
   */
  /* make (x,y) inside (-1,1) interval */
  /*x = (2 * s->x / max_x) - 1.0;
  y = (2 * s->y / max_y) - 1.0;*/
  /* make a circle of radius 0.5, centered at (0,0)->(max_x/2, max_y/2) */
  /*if( x*x + y*y < 0.25 ) g = 0.1;
  else g = 1.0;*/
  
  if(dsquare > rsquare) return(FALSE);

  d = sqrt(dsquare);

  dx /= d;
  dy /= d;
  
  f = fabs(1.0 - (d/nradius));
  
  /* include adhesion effects. I want to try adhesion
     values bigger than 1.0 */
  /* adh = 1.0 - adhes[s->type.cell_type][ttype]; */
  if ( adhes[s->ctype][ttype] > 1.0 )
      adh = adhes[s->ctype][ttype];
  else adh = 1.0 - adhes[s->ctype][ttype];
  
  *dxi = f * adh * dx * nradius;
  *dyi = f * adh * dy * nradius;

  tmpx = *dxi;
  tmpy = *dyi; 
  
  /* compute new position */
  *dxi = AniX + AniCommon * tmpx; 	/* include anisotropy effect */
  *dyi = AniY + AniCommon * tmpy; 
  
  
  /* Change for anisotropy */
  //  tmpx = *dxi;
  //  tmpy = *dyi;
  /* *dxi = *dxi * Ox - *dyi * Oy;
   *dyi = *dxi * Oy + *dyi * Ox; */

  return(TRUE);

}


/*
 * -----------------------------------------------------------------------
 * Different way to compute the forces. The neighbors here are defined
 * in a Voronoi sense
 * ----------------------------------------------------------------------- 
 */
void compute_forces_v(CELL *s, double *sumdx, double *sumdy)
{
  CELL *n, *foo;
  int numberOfNeigh;
  int i, j, k;
  double fx=0.0,fy=0.0;
	
  *sumdx = *sumdy = 0.0;
  for(i=0; i < npoints; i++){
    foo = myFacesList[i].site;
    if (foo == s){
      numberOfNeigh = myFacesList[i].faceSize;
      for(j=0; j < numberOfNeigh; j++){
	k = myFacesList[i].vorPts[j];
	n = myVoronoiList[k].sites[0];
	/* get the force between cells */
	force(s, n->ctype, n->x, n->y, &fx,&fy);
	*sumdx += fx;
	*sumdy += fy;
      }
      return;
    }
  }
}









